function [tdata,ydata] = rk4(odefun,tspan,y0,opt)
    t0  = tspan(1);
    h   = opt.step;
    tf  = tspan(2);
    y = y0;
    tdata=t0:h:tf;
    ydata=repmat(reshape(y0,1,[]),length(tdata),1);
    for i=1:length(tdata)-1
        t=tdata(i);
        k1 = h * odefun(t,y);
        k2 = h * odefun(t + h/2,y + k1/2);
        k3 = h * odefun(t + h/2,y + k2/2);
        k4 = h * odefun(t + h,y + k3);
        
        y = y + (k1 + 2*k2 + 2*k3 + k4)/6;
        if any(isnan([k1 k2 k3 k4 y]),"all")
            ydata=ydata(1:i,:);
            tdata=tdata(1:i);
            break
        end
        ydata(i+1,:)=y;
    end
end